/**********************************************
Results: HD-Mixed Sections vs LD-Mixed Sections
***********************************************/

* Clear workspace and set settings
clear all
set more off
set scheme plotplainblind


*-------------------------------------------------------------*
* Merge Line-section-level data with section-level aggregates *
*-------------------------------------------------------------*

* Load Data
use "$Data/Final/linesectionlevel_data.dta"

* Merging with section-level characteristics data
merge m:1 Line_R team new_section using "$Data/Original/Section_Aggregates.dta"
 
keep if _merge == 3 /* Egg Room and Maida (Flour) Room do not merge because no production ratings for them */
drop _merge

sort Line section Date
save "$Data/Final/linesectionlevel_datamerged.dta", replace 

*-----------------------------*
*          Switches           *
*-----------------------------*
local main_regressions     1 /* Main Regressions */
local dep_robustness       1 /* Robustness to alternative dependency thresholds */
local hindu_exposure       1 /* Heterogeneity by past exposure of Hindus */
local bjp_nrc_hindu        1 /* Heterogeneity by support for BJP/NRC among Hindus */
local shift_het            1 /* Shift heterogeneity (HD and LD pooled) */
local ritest               1 /* Permutation test */
local key_controls         1 /* Key interacted controls */
local line_clustering      1 /* Line clustering */
local baselineexp_index    1 /* Heterogeneity by baseline exposure outside work */
local spillover            1 /* Robustness to spillovers from upstream to downstream sections */

*--------------------------------------*
*         Main Regressions             *
*--------------------------------------*
if `main_regressions' == 1 {

    * Load Cleaned Data
    use "$Data/Final/linesectionlevel_datamerged.dta", clear
	
	
	*-----------------------------*
    *     Main Regressions        *
    *        Table 3            *
    *-----------------------------*

    areg rating mixed  i.line_section i.Shift, absorb(doy) cluster(line_section_team)
    eststo A1
    local obs1 = e(N)
    summarize rating if e(sample)
    local meanA = string(round(r(mean), .01))

    areg rating c.mixed#i.dep  i.line_section i.Shift, absorb(doy) cluster(line_section_team)
    eststo A2
    local obs2 = e(N)
    test 0.dep#c.mixed = 1.dep#c.mixed
    local p1 = string(round(r(p), .01))

    areg rating_dum mixed  i.line_section i.Shift, absorb(doy) cluster(line_section_team)
    eststo A3
    local obs3 = e(N)
    summarize rating_dum if e(sample)
    local meanB = string(round(r(mean), .01))

    areg rating_dum c.mixed#i.dep i.line_section i.Shift, absorb(doy) cluster(line_section_team)
    eststo A4
    local obs4 = e(N)
    test 0.dep#c.mixed = 1.dep#c.mixed
    local p2 = string(round(r(p), .01))

    * Output Table for Table 4
    estout A1 A2 A3 A4 using "$Output/Tables/tables_rating_pooledmain", style(tex) replace ///
        keep(mixed 0.dependency#c.mixed 1.dependency#c.mixed) ///
        cells(b(star fmt(%9.4f)) se(par)) ///
        nolabel collabels(none) mlabels(none) ///
        starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels( ///
            mixed                     "Mixed" ///
            0.dependency#c.mixed      "Mixed X LD" ///
            1.dependency#c.mixed      "Mixed X HD")

    * Notes for the table
    local tex " \\ \hline"
    local tex "`tex' Observations & `obs1' & `obs2' & `obs3' & `obs4'  \\"
    local tex "`tex' p(Mixed X HD = Mixed X LD) &  & `p1' &  & `p2'  \\"
    local tex "`tex' Mean Dep. Var & `meanA' & `meanA' & `meanB' & `meanB'  \\"
    local tex "`tex' Education and Tenure Controls & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Day Effects & Yes & Yes & Yes & Yes  \\"
    local tex "`tex' Shift Effects & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Line*Section Effects & Yes & Yes & Yes & Yes  \\"
    local tex "`tex' \multicolumn{5}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
    local tex " Mixed is a dummy variable coded 1 if the line-section-level team is religiously mixed."
    local tex " *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"

    esttab A1 A2 A3 A4 using "$Output/Tables/tables_rating_pooledmaina", style(tex) replace booktabs ///
        d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)" "(3)" "(4)") ///
        mgroups("Ratings" "Ratings > Median", pattern(1 0 1 0) ///
            prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
			
	
	*-----------------------------*
    *     Main Regressions        *
    *        Supervisor FE        *
	*        Table B.7            *
    *-----------------------------*

    areg rating mixed  i.line_section i.Shift i.shift_incharge, absorb(doy) cluster(line_section_team)
    eststo A1
    local obs1 = e(N)
    summarize rating if e(sample)
    local meanA = string(round(r(mean), .01))

    areg rating c.mixed#i.dep  i.line_section i.Shift i.shift_incharge, absorb(doy)  cluster(line_section_team)
    eststo A2
    local obs2 = e(N)
    test 0.dep#c.mixed = 1.dep#c.mixed
    local p1 = string(round(r(p), .01))

    areg rating_dum mixed  i.line_section i.Shift i.shift_incharge, absorb(doy)  cluster(line_section_team)
    eststo A3
    local obs3 = e(N)
    summarize rating_dum if e(sample)
    local meanB = string(round(r(mean), .01))

    areg rating_dum c.mixed#i.dep i.line_section i.Shift i.shift_incharge, absorb(doy)  cluster(line_section_team)
    eststo A4
    local obs4 = e(N)
    test 0.dep#c.mixed = 1.dep#c.mixed
    local p2 = string(round(r(p), .01))

    * Output Table for Table 4
    estout A1 A2 A3 A4 using "$Output/Tables/tables_rating_pooledmain_B7", style(tex) replace ///
        keep(mixed 0.dependency#c.mixed 1.dependency#c.mixed) ///
        cells(b(star fmt(%9.4f)) se(par)) ///
        nolabel collabels(none) mlabels(none) ///
        starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels( ///
            mixed                     "Mixed" ///
            0.dependency#c.mixed      "Mixed X LD" ///
            1.dependency#c.mixed      "Mixed X HD")

    * Notes for the table
    local tex " \\ \hline"
    local tex "`tex' Observations & `obs1' & `obs2' & `obs3' & `obs4'  \\"
    local tex "`tex' p(Mixed X HD = Mixed X LD) &  & `p1' &  & `p2'  \\"
    local tex "`tex' Mean Dep. Var & `meanA' & `meanA' & `meanB' & `meanB'  \\"
    local tex "`tex' Education and Tenure Controls & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Day Effects & Yes & Yes & Yes & Yes  \\"
    local tex "`tex' Shift Effects & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Line*Section Effects & Yes & Yes & Yes & Yes  \\"
	local tex "`tex' Supervisor Effects & Yes & Yes & Yes & Yes  \\"
    local tex "`tex' \multicolumn{5}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
    local tex " Mixed is a dummy variable coded 1 if the line-section-level team is religiously mixed."
    local tex " *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"

    esttab A1 A2 A3 A4 using "$Output/Tables/tables_rating_pooledmaina_B7", style(tex) replace booktabs ///
        d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)" "(3)" "(4)") ///
        mgroups("Ratings" "Ratings > Median", pattern(1 0 1 0) ///
            prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
			
	


    *--------------------------------*
    *     Main Regressions           *
	*  Schooling and Tenure controls *
    *        Table B.8               *
    *--------------------------------*

    areg rating mixed school_years tenure i.line_section i.Shift, absorb(doy) cluster(line_section_team)
    eststo A1
    local obs1 = e(N)
    summarize rating if e(sample)
    local meanA = string(round(r(mean), .01))

    areg rating c.mixed#i.dep school_years tenure i.line_section i.Shift, absorb(doy) cluster(line_section_team)
    eststo A2
    local obs2 = e(N)
    test 0.dep#c.mixed = 1.dep#c.mixed
    local p1 = string(round(r(p), .01))

    areg rating_dum mixed school_years tenure i.line_section i.Shift, absorb(doy) cluster(line_section_team)
    eststo A3
    local obs3 = e(N)
    summarize rating_dum if e(sample)
    local meanB = string(round(r(mean), .01))

    areg rating_dum c.mixed#i.dep school_years tenure i.line_section i.Shift, absorb(doy) cluster(line_section_team)
    eststo A4
    local obs4 = e(N)
    test 0.dep#c.mixed = 1.dep#c.mixed
    local p2 = string(round(r(p), .01))

    * Output Table for Table B.6
    estout A1 A2 A3 A4 using "$Output/Tables/tables_rating_pooled_B8", style(tex) replace ///
        keep(mixed 0.dependency#c.mixed 1.dependency#c.mixed) ///
        cells(b(star fmt(%9.4f)) se(par)) ///
        nolabel collabels(none) mlabels(none) ///
        starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels( ///
            mixed                     "Mixed" ///
            0.dependency#c.mixed      "Mixed X LD" ///
            1.dependency#c.mixed      "Mixed X HD")

    * Notes for the table
    local tex " \\ \hline"
    local tex "`tex' Observations & `obs1' & `obs2' & `obs3' & `obs4'  \\"
    local tex "`tex' p(Mixed X HD = Mixed X LD) &  & `p1' &  & `p2'  \\"
    local tex "`tex' Mean Dep. Var & `meanA' & `meanA' & `meanB' & `meanB'  \\"
    local tex "`tex' Education and Tenure Controls & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Day Effects & Yes & Yes & Yes & Yes  \\"
    local tex "`tex' Shift Effects & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Line*Section Effects & Yes & Yes & Yes & Yes  \\"
    local tex "`tex' \multicolumn{5}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
    local tex " Mixed is a dummy variable coded 1 if the line-section-level team is religiously mixed."
    local tex " *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"

    esttab A1 A2 A3 A4 using "$Output/Tables/tables_rating_pooleda_B8", style(tex) replace booktabs ///
        d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)" "(3)" "(4)") ///
        mgroups("Ratings" "Ratings > Median", pattern(1 0 1 0) ///
            prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))

    *-----------------------------*
    *    HD and LD Separately     *
    *       Table B.5             *
    *-----------------------------*

    * HD Sections
    areg rating mixed i.line_section i.Shift if dep == 1, absorb(doy) cluster(line_section_team)
    eststo B1
    local obs1 = e(N)
    summarize rating if e(sample)
    local meanA = string(round(r(mean), .01))

    areg rating_dum mixed i.line_section i.Shift if dep == 1, absorb(doy) cluster(line_section_team)
    eststo B2
    local obs2 = e(N)
    summarize rating_dum if e(sample)
    local meanB = string(round(r(mean), .01))

    * LD Sections
    areg rating mixed school_years tenure i.line_section i.Shift if dep == 0, absorb(doy) cluster(line_section_team)
    eststo B3
    local obs3 = e(N)
    summarize rating if e(sample)
    local meanC = string(round(r(mean), .01))

    areg rating_dum mixed school_years tenure i.line_section i.Shift if dep == 0, absorb(doy) cluster(line_section_team)
    eststo B4
    local obs4 = e(N)
    summarize rating_dum if e(sample)
    local meanD = string(round(r(mean), .01))

    * Output Table for Table B.5
    estout B1 B2 B3 B4 using "$Output/Tables/tables_rating_HDvsLD", style(tex) replace ///
        keep(mixed) ///
        cells(b(star fmt(%9.4f)) se(par)) ///
        nolabel collabels(none) mlabels(none) ///
        starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels(mixed "Mixed")

    * Notes for the table
    local tex " \\ \hline"
    local tex "`tex' Observations & `obs1' & `obs2' & `obs3' & `obs4'  \\"
    local tex "`tex' Mean Dep. Var & `meanA' & `meanB' & `meanC' & `meanD'  \\"
    local tex "`tex' Education and Tenure Controls & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Day Effects & Yes & Yes & Yes & Yes  \\"
    local tex "`tex' Shift Effects & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Line*Section Effects & Yes & Yes & Yes & Yes  \\"
    local tex "`tex' \multicolumn{5}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
    local tex " Mixed is a dummy variable coded 1 if the line-section-level team is religiously mixed."
    local tex " *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"

    esttab B1 B2 B3 B4 using "$Output/Tables/tables_rating_HDvsLDa", style(tex) replace booktabs ///
        d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)" "(3)" "(4)") ///
        mgroups("Ratings" "Ratings > Median", pattern(1 1 1 1) ///
            prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))

*--------------------------------------*
*   Effects During Days of Violent     *
*   Rioting in West Bengal and Delhi   *
*   over the Citizenship Amendment Act *
*--------------------------------------*


*--------------------------------------*
*     Create Violence Indicator        *
*--------------------------------------*

* Generate 'violence' dummy variable
generate violence = 0
replace violence = 1 if Date >= date("12dec2019", "DMY") & Date <= date("18dec2019", "DMY")
* CAA announced and subsequent violence in Howrah
replace violence = 1 if Date >= date("23feb2020", "DMY") & Date <= date("29feb2020", "DMY")
* CAA related riots in New Delhi

*--------------------------------------*
*       High Dependency (HD)           *
*             Sections                 *
*--------------------------------------*

* Model 1: Without interaction
reghdfe rating mixed if dep == 1, absorb(doy i.line_section i.Shift) cluster(line_section_team)
estimates store A1
local obs1 = e(N)
summarize rating if e(sample)
local meanA = string(round(r(mean), .01))

* Model 2: With interaction between mixed and violence
reghdfe rating c.mixed#i.violence if dep == 1, absorb(doy i.line_section i.Shift) cluster(line_section_team)
estimates store A2
local obs2 = e(N)
summarize rating if e(sample)
local meanB = string(round(r(mean), .01))

* Test for difference in mixed effect between violence periods
test 0.violence#c.mixed = 1.violence#c.mixed

*--------------------------------------*
*       Low Dependency (LD)            *
*             Sections                 *
*--------------------------------------*

* Model 3: Without interaction
reghdfe rating mixed if dep == 0, absorb(doy i.line_section i.Shift) cluster(line_section_team) 
estimates store A3
local obs3 = e(N)
summarize rating if e(sample)
local meanC = string(round(r(mean), .01))

* Model 4: With interaction between mixed and violence
reghdfe rating c.mixed#i.violence if dep == 0, absorb(doy i.line_section i.Shift) cluster(line_section_team)
estimates store A4
local obs4 = e(N)
summarize rating if e(sample)
local meanD = string(round(r(mean), .01))

* Test for difference in mixed effect between violence periods
test 0.violence#c.mixed = 1.violence#c.mixed

*--------------------------------------*
*             Output Table             *
*--------------------------------------*

estout A1 A2 A3 A4 using  "$Output/Tables/tables_rating_RV", style(tex) replace ///
	keep(mixed 0.violence#c.mixed 1.violence#c.mixed) ///
	cells(b(star fmt(%9.4f)) se(par)) ///
	nolabel collabels(none) mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) ///
	varlabels(mixed  "Mixed" ///
			 0.violence#c.mixed  "Mixed X No Violence" ///
			 1.violence#c.mixed  "Mixed X Violence")
			 
			 
local tex " \\ \hline"
local tex "`tex' Observations & `obs1' & `obs2' & `obs3' & `obs4'   \\"
local tex "`tex' Mean Dep. Var & `meanA' & `meanB' & `meanC' & `meanD'  \\"

local tex "`tex' Education and Tenure Controls & Yes & Yes & Yes & Yes \\"
local tex "`tex' Day Effects & Yes & Yes & Yes & Yes \\"
local tex "`tex' Shift Effects & Yes & Yes & Yes & Yes \\"
local tex "`tex' Line*Section Effects & Yes & Yes & Yes & Yes \\"
local tex "`tex' \multicolumn{5}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
local tex "`tex' Mixed  is  a  dummy  variable  coded  1  if  the  line-section-level  team  is  religiously mixed." 
local tex "`tex' *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"
	
esttab A1 A2 A3 A4 using "$Output/Tables/tables_rating_RVa", style(tex) replace booktabs ///
	d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
	mtitles("(1)" "(2)" "(3)" "(4)") ///
	mgroups("HD Sections" "LD Sections", pattern(1 0 1 0) ///
	prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
	
     
*--------------------------------------*
*     Convergence Over Time Analysis   *
*--------------------------------------*


*--------------------------------------*
*          Create Day Bins             *
*--------------------------------------*

* Create day bins to show convergence over time
egen day_bins = cut(int_days), at(0(25)125)
replace day_bins=100 if mi(day_bins)
*--------------------------------------*
*       High Dependency (HD)           *
*             Sections                 *
*--------------------------------------*

* Model A1: Without interaction
reghdfe rating mixed school_years tenure if dep == 1, absorb(doy i.line_section i.Shift) cluster(line_section_team) 
estimates store A1
local obs1 = e(N)
summarize rating if e(sample)
local meanA = string(round(r(mean), .01))

* Model A2: With interaction between mixed and day_bins
reghdfe rating c.mixed#i.day_bins school_years tenure if dep == 1, absorb(doy i.line_section i.Shift) cluster(line_section_team)
estimates store A2
local obs2 = e(N)

*--------------------------------------*
*       Low Dependency (LD)            *
*             Sections                 *
*--------------------------------------*

* Model A3: Without interaction
reghdfe rating mixed school_years tenure if dep == 0, absorb(doy i.line_section i.Shift) cluster(line_section_team) 
estimates store A3
local obs3 = e(N)
summarize rating if e(sample)
local meanB = string(round(r(mean), .01))

* Model A4: With interaction between mixed and day_bins
reghdfe rating c.mixed#i.day_bins school_years tenure if dep == 0, absorb(doy i.line_section i.Shift) cluster(line_section_team) 
estimates store A4
local obs4 = e(N)

*--------------------------------------*
*             Output Table             *
*--------------------------------------*

estout A1 A2 A3 A4 using  "$Output/Tables/tables_rating_att", style(tex) replace ///
	keep(mixed 0.day_bins#c.mixed  25.day_bins#c.mixed 50.day_bins#c.mixed 75.day_bins#c.mixed 100.day_bins#c.mixed) ///
	cells(b(star fmt(%9.4f)) se(par)) ///
	nolabel collabels(none) mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) ///
	varlabels(mixed "Mixed" ///
	             0.day_bins#c.mixed  "Mixed X 0-25 days" ///
				 25.day_bins#c.mixed  "Mixed X 26-50 days" ///
				50.day_bins#c.mixed "Mixed X 51-75 days" ///
				75.day_bins#c.mixed "Mixed X 76-100 days" ///
				100.day_bins#c.mixed "Mixed X > 100 days") 
				
				
				
local tex " \\ \hline"
local tex "`tex' Observations & `obs1' & `obs2' & `obs3' & `obs4'  \\"
local tex "`tex' Mean Dep. Var & `meanA' & `meanA' & `meanB'' & `meanB'  \\"
local tex "`tex' Education and Tenure Controls & Yes & Yes & Yes & Yes \\"
local tex "`tex' Day Effects & Yes & Yes & Yes & Yes  \\"
local tex "`tex' Shift Effects & Yes & Yes & Yes & Yes \\"
local tex "`tex' Line*Section Effects & Yes & Yes & Yes & Yes  \\"
local tex "`tex' \multicolumn{5}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
local tex "`tex' Mixed  is  a  dummy  variable  coded  1  if  the  line-section-level  team  is  religiously mixed." 
local tex "`tex' *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"
	
esttab A1 A2 A3 A4 using "$Output/Tables/tables_rating_atta", style(tex) replace booktabs ///
	d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
	mtitles("(1)" "(2)" "(3)" "(4)") ///
	mgroups("Ratings", pattern(1 0 0 0) ///
	prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
	




}

*--------------------------------------*
*         Key Controls (Table B.9)     *
*--------------------------------------*
if `key_controls' == 1 {

    * Load Cleaned Data
    use "$Data/Final/linesectionlevel_datamerged.dta", clear

    
    * Generate interacted controls
    generate mixedXcountw_std        = mixed * count_workers_std
    generate LDXtenure_s_std         = LD * tenure_skilled_std
    generate HDXtenure_s_std         = HD * tenure_skilled_std
    generate LDXtenure_us_std        = LD * tenure_unskilled_std
    generate HDXtenure_us_std        = HD * tenure_unskilled_std
    generate LDXschool_std           = LD * school_years_std
    generate HDXschool_std           = HD * school_years_std
    generate mixedXhs_hindu_std      = mixed * share_hskilled_hindus_std
    generate mixedXhs_muslim_std     = mixed * share_hskilled_muslims_std
    generate mixedXtenure_s_std      = mixed * tenure_skilled_std
    generate mixedXtenure_us_std     = mixed * tenure_unskilled_std
	generate mixedXage_std           = mixed * age_std
	
	* Save data for starbility plot
	save $Data/Final/linesectionlevel_starbility.dta, replace

    *-----------------------------*
    *   Main Regressions with     *
    *     Key Controls            *
    *        Table 6              *
    *-----------------------------*

    * Regression 1
    reghdfe rating mixedXLD mixedXHD, absorb(doy line_section Shift) cluster(line_section_team)
    eststo A1
    local obs1 = e(N)
    summarize rating if e(sample)
    local mean1 = string(round(r(mean), .01))
    test mixedXHD = mixedXLD

    * Regression 2: Add mixedXcountw_std
    reghdfe rating mixedXLD mixedXHD mixedXcountw_std, absorb(doy line_section Shift) cluster(line_section_team)
    eststo A2
    local obs2 = e(N)
    summarize rating if e(sample)
    local mean2 = string(round(r(mean), .01))
    test mixedXHD = mixedXLD

    * Regression 3: Add LDXschool_std and HDXschool_std
    reghdfe rating mixedXLD mixedXHD mixedXcountw_std LDXschool_std HDXschool_std, absorb(doy line_section Shift) cluster(line_section_team)
    eststo A3
    local obs3 = e(N)
    summarize rating if e(sample)
    local mean3 = string(round(r(mean), .01))
    test mixedXHD = mixedXLD

    * Regression 4: Add tenure controls
    reghdfe rating mixedXLD mixedXHD mixedXcountw_std LDXschool_std HDXschool_std ///
        LDXtenure_s_std HDXtenure_s_std LDXtenure_us_std HDXtenure_us_std ///
        mixedXtenure_us_std mixedXtenure_s_std, absorb(doy line_section Shift) cluster(line_section_team)
    eststo A4
    local obs4 = e(N)
    summarize rating if e(sample)
    local mean4 = string(round(r(mean), .01))
    test mixedXHD = mixedXLD

    * Regression 5: Add high_skilled variables
    reghdfe rating mixedXLD mixedXHD mixedXcountw_std LDXtenure_s_std HDXtenure_s_std ///
        LDXtenure_us_std HDXtenure_us_std mixedXtenure_us_std mixedXtenure_s_std ///
        LDXschool_std HDXschool_std mixedXhs_hindu_std mixedXhs_muslim_std ///
       age_std mixedXage_std, absorb(doy line_section Shift) cluster(line_section_team)
    eststo A5
    local obs5 = e(N)
    summarize rating if e(sample)
    local mean5 = string(round(r(mean), .01))
    test mixedXHD = mixedXLD

    * Output Table for Table B.9
    estout A1 A2 A3 A4 A5 using "$Output/Tables/tables_rating_pooledcon", style(tex) replace ///
        keep(mixedXLD mixedXHD mixedXcountw_std LDXtenure_s_std HDXtenure_s_std ///
            LDXtenure_us_std HDXtenure_us_std LDXschool_std HDXschool_std ///
            mixedXhs_hindu_std mixedXhs_muslim_std mixedXtenure_us_std mixedXtenure_s_std) ///
        cells(b(star fmt(%9.4f)) se(par)) ///
        nolabel collabels(none) mlabels(none) ///
        starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels( ///
            mixedXLD                "Mixed $\times$ LD" ///
            mixedXHD                "Mixed $\times$ HD" ///
            mixedXcountw_std        "Mixed $\times$ Group Size" ///
            LDXtenure_s_std         "LD $\times$ Tenure (Skilled)" ///
            HDXtenure_s_std         "HD $\times$ Tenure (Skilled)" ///
            LDXtenure_us_std        "LD $\times$ Tenure (Unskilled)" ///
            HDXtenure_us_std        "HD $\times$ Tenure (Unskilled)" ///
            mixedXtenure_us_std     "Mixed $\times$ Tenure (Unskilled)" ///
            mixedXtenure_s_std      "Mixed $\times$ Tenure (Skilled)" ///
            LDXschool_std           "LD $\times$ Schooling" ///
            HDXschool_std           "HD $\times$ Schooling" ///
            mixedXhs_hindu_std      "Mixed $\times$ Share Skilled (Hindu)" ///
            mixedXhs_muslim_std     "Mixed $\times$ Share Skilled (Muslim)")

    * Notes for the table
    local tex " \\ \hline"
    local tex "`tex' Day F.E. & Yes & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Shift F.E. & Yes & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Line $\times$ Section F.E. & Yes & Yes & Yes & Yes & Yes \\"
    local tex "`tex' Mean Dep Var. & `mean1' & `mean2' & `mean3' & `mean4'  & `mean5'  \\"
    local tex "`tex' N & `obs1' & `obs2' & `obs3' & `obs4' & `obs5' \\"
    local tex "`tex' \multicolumn{6}{p{11cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
    local tex " Mixed is a dummy variable coded 1 if the line-section-level team is religiously mixed."
    local tex " *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"

    esttab A1 A2 A3 A4 A5 using "$Output/Tables/tables_rating_pooledcona", style(tex) replace booktabs ///
        d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)" "(3)" "(4)" "(5)") ///
        mgroups("Ratings", pattern(1 1 1 1 1) ///
            prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))

}


*--------------------------------------*
*              Permuation Test         *
*--------------------------------------*
if `ritest' == 1{

*--------------------------------------*
*              Table B18               *
*--------------------------------------*

	
* Load Cleaned Data
use "$Data/Final/linesectionlevel_datamerged.dta", clear

set seed 3585051

* Creating a single variable for treatment for randomization inference.
gen treatment = 0
replace treatment = 1 if mixedXLD == 1
replace treatment = 2 if mixedXHD == 1

* Mixed regression
{
	* Regress with rating
		reghdfe rating i.mixed, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		* Store the regression result
		estimates store modelA

		* Extract coefficients and standard errors, then round to three decimal places
		local mixed_coefA = string(round(e(b)[1,"1.mixed"], 0.001), "%9.3f")
		local obsA = e(N)
		local r2A = string(round(e(r2_a), 0.001), "%9.3f")

		su rating if e(sample) == 1
		local meanA = string(round(r(mean), .01))

		* Calculate p-values and append significance stars
		local mixed_tA = e(b)[1,"1.mixed"]/sqrt(e(V)[2,2])
		local mixed_pvalA = string(round(2*ttail(e(df_r),abs(`mixed_tA')), 0.001), "%9.3f")
		local mixed_starsA = cond(`mixed_pvalA' < 0.01, "***", cond(`mixed_pvalA' < 0.05, "**", cond(`mixed_pvalA' < 0.10, "*", "")))
		
		* Combine coefficients and stars
		local mixed_resultA = "`mixed_coefA'`mixed_starsA'"
		local mixed_resultA = "`mixed_coefA'`mixed_starsA'"

		* Perform ritest for mixed and store p-values
		ritest mixed _b[1.mixed], cluster(line_section_team) strata(line_section) reps(2000) nodots: reghdfe rating i.mixed, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		local mixed_riA = string(round(el(r(p), 1, 1), 0.001), "%9.3f")
	
	* Regress with rating dummy
		reghdfe rating_dum i.mixed, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		* Store the regression result
		estimates store modelB

		* Extract coefficients and standard errors, then round to three decimal places
		local mixed_coefB = string(round(e(b)[1,"1.mixed"], 0.001), "%9.3f")
		local obsB = e(N)
		local r2B = string(round(e(r2_a), 0.001), "%9.3f")

		su rating_dum if e(sample) == 1
		local meanB = string(round(r(mean), .01))

		* Calculate p-values and append significance stars
		local mixed_tB = e(b)[1,"1.mixed"]/sqrt(e(V)[2,2])
		local mixed_pvalB = string(round(2*ttail(e(df_r),abs(`mixed_tB')), 0.001), "%9.3f")
		local mixed_starsB = cond(`mixed_pvalB' < 0.01, "***", cond(`mixed_pvalB' < 0.05, "**", cond(`mixed_pvalB' < 0.10, "*", "")))
		
		* Combine coefficients and stars
		local mixed_resultB = "`mixed_coefB'`mixed_starsB'"
		local mixed_resultB = "`mixed_coefB'`mixed_starsB'"

		* Perform ritest for mixed and store p-values
		ritest mixed _b[1.mixed], cluster(line_section_team) strata(line_section) reps(2000) nodots: reghdfe rating_dum i.mixed, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		local mixed_riB = string(round(el(r(p), 1, 1), 0.001), "%9.3f")
}


* Mixed interaction regression
{
	* Regress with rating
		reghdfe rating i.treatment, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		* Store the regression result
		estimates store modelXA

		* Extract coefficients and standard errors, then round to three decimal places
		local mixedXLD_coefA = string(round(e(b)[1,"1.treatment"], 0.001), "%9.3f")
		local mixedXHD_coefA = string(round(e(b)[1,"2.treatment"], 0.001), "%9.3f")
		local obsXA = e(N)
		local r2XA = string(round(e(r2_a), 0.001), "%9.3f")

		su rating if e(sample) == 1
		local meanXA = string(round(r(mean), .01))

		* Calculate p-values and append significance stars
		local mixedXLD_tA = e(b)[1,"1.treatment"]/sqrt(e(V)[2,2])
		local mixedXLD_pvalA = string(round(2*ttail(e(df_r),abs(`mixedXLD_tA')), 0.001), "%9.3f")
		local mixedXLD_starsA = cond(`mixedXLD_pvalA' < 0.01, "***", cond(`mixedXLD_pvalA' < 0.05, "**", cond(`mixedXLD_pvalA' < 0.10, "*", "")))
		local mixedXHD_tA = e(b)[1,"2.treatment"]/sqrt(e(V)[3,3])
		local mixedXHD_pvalA = string(round(2*ttail(e(df_r),abs(`mixedXHD_tA')), 0.001), "%9.3f")
		local mixedXHD_starsA = cond(`mixedXHD_pvalA' < 0.01, "***", cond(`mixedXHD_pvalA' < 0.05, "**", cond(`mixedXHD_pvalA' < 0.10, "*", "")))

		* Combine coefficients and stars
		local mixedXLD_resultA = "`mixedXLD_coefA'`mixedXLD_starsA'"
		local mixedXHD_resultA = "`mixedXHD_coefA'`mixedXHD_starsA'"

		* Perform ritest for mixedXLD and store p-values
		ritest treatment _b[1.treatment], fixlevels(2) cluster(line_section_team) strata(line_section) reps(2000) nodots: reghdfe rating i.treatment, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		local mixedXLD_riA = string(round(el(r(p), 1, 1), 0.001), "%9.3f")

		* Perform ritest for mixedXHD and store p-values
		ritest treatment _b[2.treatment], fixlevels(1) cluster(line_section_team) strata(line_section) reps(2000) nodots: reghdfe rating i.treatment, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		local mixedXHD_riA = string(round(el(r(p), 1, 1), 0.001), "%9.3f")
	
	* Regress with rating dummy
		reghdfe rating_dum i.treatment, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		* Store the regression result
		estimates store modelXB

		* Extract coefficients and standard errors, then round to three decimal places
		local mixedXLD_coefB = string(round(e(b)[1,"1.treatment"], 0.001), "%9.3f")
		local mixedXHD_coefB = string(round(e(b)[1,"2.treatment"], 0.001), "%9.3f")
		local obsXB = e(N)
		local r2XB = string(round(e(r2_a), 0.001), "%9.3f")

		su rating_dum if e(sample) == 1
		local meanXB = string(round(r(mean), .01))

		* Calculate p-values and append significance stars
		local mixedXLD_tB = e(b)[1,"1.treatment"]/sqrt(e(V)[2,2])
		local mixedXLD_pvalB = string(round(2*ttail(e(df_r),abs(`mixedXLD_tB')), 0.001), "%9.3f")
		local mixedXLD_starsB = cond(`mixedXLD_pvalB' < 0.01, "***", cond(`mixedXLD_pvalB' < 0.05, "**", cond(`mixedXLD_pvalB' < 0.10, "*", "")))
		local mixedXHD_tB = e(b)[1,"2.treatment"]/sqrt(e(V)[3,3])
		local mixedXHD_pvalB = string(round(2*ttail(e(df_r),abs(`mixedXHD_tB')), 0.001), "%9.3f")
		local mixedXHD_starsB = cond(`mixedXHD_pvalB' < 0.01, "***", cond(`mixedXHD_pvalB' < 0.05, "**", cond(`mixedXHD_pvalB' < 0.10, "*", "")))

		* Combine coefficients and stars
		local mixedXLD_resultB = "`mixedXLD_coefB'`mixedXLD_starsB'"
		local mixedXHD_resultB = "`mixedXHD_coefB'`mixedXHD_starsB'"

		* Perform ritest for mixedXLD and store p-values
		ritest treatment _b[1.treatment], fixlevels(2) cluster(line_section_team) strata(line_section) reps(2000) nodots: reghdfe rating_dum i.treatment, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		local mixedXLD_riB = string(round(el(r(p), 1, 1), 0.001), "%9.3f")

		* Perform ritest for mixedHD and store p-values
		ritest treatment _b[2.treatment], fixlevels(1) cluster(line_section_team) strata(line_section) reps(2000) nodots: reghdfe rating_dum i.treatment, absorb(doy i.line_section i.Shift) cluster(line_section_team)

		local mixedXHD_riB = string(round(el(r(p), 1, 1), 0.001), "%9.3f")
}

* Open a .tex file to store the results
file open results using $Output/Tables/tables_ritest_rating.tex, write replace

* Write the results to the .tex file
file write results "Mixed & `mixed_resultA' &	& `mixed_resultB' &  \\" _n
file write results "  & (`mixed_pvalA') &	& (`mixed_pvalB') &   \\" _n
file write results "  & [`mixed_riA'] &	& [`mixed_riB'] &   \\" _n
file write results "Mixed $\times$ LD &	& `mixedXLD_resultA' &	& `mixedXLD_resultB' \\" _n
file write results "  & & (`mixedXLD_pvalA') & & (`mixedXLD_pvalB') \\" _n
file write results "  & & [`mixedXLD_riA'] & & [`mixedXLD_riB'] \\" _n
file write results "Mixed $\times$ HD &	& `mixedXHD_resultA' &	& `mixedXHD_resultB' \\" _n
file write results "  & & (`mixedXHD_pvalA') & & (`mixedXHD_pvalB') \\" _n
file write results "  & & [`mixedXHD_riA'] & & [`mixedXHD_riB'] \\" _n
file close results


local tex " \\ \hline"
local tex "`tex' Day F.E. & Yes & Yes & Yes & Yes \\"
local tex "`tex' Shift F.E. & Yes & Yes & Yes & Yes \\"
local tex "`tex' Line $\times$ Section F.E. & Yes & Yes & Yes & Yes \\"
local tex "`tex' Mean Dep Var. & `meanA' & `meanXA' & `meanB' & `meanXB' \\"
local tex "`tex' N & `obsA' & `obsXA' & `obsB' & `obsXB' \\"
local tex "`tex' Adj. $ R^2$ & `r2A' & `r2XA' & `r2B' & `r2XB' \\"
local tex "`tex' \multicolumn{5}{l}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
local tex "`tex' P-values are reported in parentheses."
local tex "`tex' Mixed  is  a  dummy  variable  coded  1 if  the  line-section-level  team  is  religiously mixed." 
local tex "`tex' *** p$<$0.01, ** p$<$0.05, * p$<$0.1.} \\ \end{tabular} }"
	
esttab modelA modelXA modelB modelXB ///
	using "$Output/Tables/tables_ritest_ratinga", style(tex) replace booktabs ///
	d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
	mtitles("(1)" "(2)" "(3)" "(4)") ///
	mgroups("Ratings", pattern(1) ///
	prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
}


*--------------------------------------*
*       Line Clustering Analysis       *
*--------------------------------------*
if `line_clustering' == 1{

*------------*
* Table B19  *
*------------*

   * Load Cleaned Data
    use "$Data/Final/linesectionlevel_datamerged.dta", clear


* Regressions
 areg rating mixed  i.line_section i.Shift, absorb(doy) cluster(line_team)
 eststo A1
 local A1_pval = string(round(r(table)[4,1], 0.001), "%9.3f")
 local A1_coef = string(round(r(table)[1,1], 0.001), "%9.3f")
 local A1_stars = cond(`A1_pval' < 0.01, "***", cond(`A1_pval' < 0.05, "**", cond(`A1_pval' < 0.10, "*", "")))
 local A1_result = "`A1_coef'`A1_stars'"
 
 local obs1 = e(N)
 su rating if e(sample) == 1
 local meanA = string(round(r(mean), .01))
 
 wildbootstrap areg rating mixed  i.line_section i.Shift, absorb(doy) cluster(line_team)   coefficients(mixed) rseed(2315)
 local A1_pval_wbs = string(round(r(table)[3,1], 0.001), "%9.3f")
 
 areg rating mixedXLD mixedXHD i.line_section i.Shift , absorb(doy) cluster(line_team)
 matrix list r(table)
 eststo A2
 local LD_A2_pval = string(round(r(table)[4,1], 0.001), "%9.3f")
 local LD_A2_coef = string(round(r(table)[1,1], 0.001), "%9.3f")
 local LD_A2_stars = cond(`LD_A2_pval' < 0.01, "***", cond(`LD_A2_pval' < 0.05, "**", cond(`LD_A2_pval' < 0.10, "*", "")))
 local LD_A2_result = "`LD_A2_coef'`LD_A2_stars'"
 
 local HD_A2_pval = string(round(r(table)[4,2], 0.001), "%9.3f")
 local HD_A2_coef = string(round(r(table)[1,2], 0.001), "%9.3f")
 local HD_A2_stars = cond(`HD_A2_pval' < 0.01, "***", cond(`HD_A2_pval' < 0.05, "**", cond(`HD_A2_pval' < 0.10, "*", "")))
 local HD_A2_result = "`HD_A2_coef'`HD_A2_stars'"
 
 local obs2 = e(N)
 test mixedXLD = mixedXHD
 local p1 = string(round(r(p), .01))
 
 wildbootstrap areg rating mixedXLD mixedXHD i.line_section i.Shift, absorb(doy) cluster( line_team)  coefficients(mixedXLD mixedXHD) rseed(11891)
 local LD_A2_pval_wbs = string(round(r(table)[3,1], 0.001), "%9.3f")
 local HD_A2_pval_wbs = string(round(r(table)[3,2], 0.001), "%9.3f")
 

 
 

 areg rating_dum mixed  i.line_section i.Shift, absorb(doy) cluster(line_team)
 eststo A3
 local A3_pval = string(round(r(table)[4,1], 0.001), "%9.3f")
 local A3_coef = string(round(r(table)[1,1], 0.001), "%9.3f")
 local A3_stars = cond(`A3_pval' < 0.01, "***", cond(`A3_pval' < 0.05, "**", cond(`A3_pval' < 0.10, "*", "")))
 local A3_result = "`A3_coef'`A3_stars'"

 local obs3 = e(N)
 su rating_dum if e(sample) == 1
 local meanB = string(round(r(mean), .01))
 
 wildbootstrap areg rating_dum mixed  i.line_section i.Shift, absorb(doy) cluster(line_team)   coefficients(mixed) rseed(1199)
  local A3_pval_wbs = string(round(r(table)[3,1], 0.001), "%9.3f")

 
 areg rating_dum mixedXLD mixedXHD i.line_section i.Shift, absorb(doy) cluster(line_team)
 eststo A4
 local LD_A4_pval = string(round(r(table)[4,1], 0.001), "%9.3f")
 local LD_A4_coef = string(round(r(table)[1,1], 0.001), "%9.3f")
 local LD_A4_stars = cond(`LD_A4_pval' < 0.01, "***", cond(`LD_A4_pval' < 0.05, "**", cond(`LD_A4_pval' < 0.10, "*", "")))
 local LD_A4_result = "`LD_A4_coef'`LD_A4_stars'"
 
 local HD_A4_pval = string(round(r(table)[4,2], 0.001), "%9.3f")
 local HD_A4_coef = string(round(r(table)[1,2], 0.001), "%9.3f")
 local HD_A4_stars = cond(`HD_A4_pval' < 0.01, "***", cond(`HD_A4_pval' < 0.05, "**", cond(`HD_A4_pval' < 0.10, "*", "")))
 local HD_A4_result = "`HD_A4_coef'`HD_A4_stars'"
 
 local obs4 = e(N)
 test mixedXLD = mixedXHD
 local p2 = string(round(r(p), .01))

 wildbootstrap areg rating_dum mixedXLD mixedXHD i.line_section i.Shift, absorb(doy) cluster(line_team)  coefficients(mixedXLD mixedXHD) rseed(11891) 
 local LD_A4_pval_wbs = string(round(r(table)[3,1], 0.001), "%9.3f")
 local HD_A4_pval_wbs = string(round(r(table)[3,2], 0.001), "%9.3f")

* Open a .tex file to store the results
file open results using $Output/Tables/tables_rating_pooled, write replace

* Write the results to the .tex file
file write results "Mixed & `A1_result' &  & `A3_result' & \\" _n
file write results "  & (`A1_pval') &  & (`A3_pval') & \\" _n
file write results "  & [`A1_pval_wbs'] &  & [`A3_pval_wbs'] &  \\" _n
file write results "Mixed $\times$ LD &  & `LD_A2_result' &  & `LD_A4_result' \\" _n
file write results "  &  & (`LD_A2_pval') &  & (`LD_A4_pval')  \\" _n
file write results "  &  & [`LD_A2_pval_wbs'] &  & [`LD_A4_pval_wbs'] \\" _n
file write results "Mixed $\times$ HD &  & `HD_A2_result' &  & `HD_A4_result'  \\" _n
file write results "  &  & (`HD_A2_pval') &  & (`HD_A4_pval')  \\" _n
file write results "  &  & [`HD_A2_pval_wbs'] &  & [`HD_A4_pval_wbs'] \\" _n
file close results

			 
local tex " \\ \hline"
local tex "`tex' p(Mixed $\times$ HD = Mixed $\times$ LD) &  & `p1' &  & `p2'  \\"
local tex "`tex' Day F.E. & Yes & Yes & Yes & Yes  \\"
local tex "`tex' Shift F.E. & Yes & Yes & Yes & Yes \\"
local tex "`tex' Line $\times$ Section F.E. & Yes & Yes & Yes & Yes  \\"
local tex "`tex' Mean Dep Var. & `meanA' & `meanA' & `meanB' & `meanB'  \\"
local tex "`tex' N & `obs1' & `obs2' & `obs3' & `obs4'  \\"
local tex "`tex' \multicolumn{5}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
local tex "`tex' Mixed  is  a  dummy  variable  coded  1  if  the  line-section-level  team  is  religiously mixed." 
local tex "`tex' *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"
	
esttab A1 A2 A3 A4 using "$Output/Tables/tables_rating_pooleda", style(tex) replace booktabs ///
	d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
	mtitles("(1)" "(2)" "(3)" "(4)") ///
	mgroups("Ratings" "Ratings > Median", pattern(1 0 1 0) ///
	prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
	
}


*--------------------------------------*
*     Hindu Exposure Analysis          *
*--------------------------------------*
if `hindu_exposure' == 1 {
    
    
    * Set memory and options
    set more off
    clear all
    set segmentsize 3g
    set scheme plotplainblind
    
    * Load Cleaned Data
    use "$Data/Final/linesectionlevel_datamerged.dta", clear
    
    ****************************************
    ** Mixing Treatment Based on Contact **
    **** that Hindus Had at Baseline *****
    ***************************************
    
    * Generate variables for high and low contact Hindus
    gen mixed_highc_hindus = mixed
    quietly summarize mean_contact_hindu_l, detail
    local p50 = r(p50)
    
    replace mixed_highc_hindus = . if mean_contact_hindu_l < `p50' & mixed == 1  // Hindus had high exposure to Muslims
    
    gen mixed_lowc_hindus = mixed
    replace mixed_lowc_hindus = . if mean_contact_hindu_l >= `p50' & mixed == 1  // Hindus had low exposure to Muslims
    
    * Create day bins
    egen day_bins = cut(int_days), at(0(20)120)
    replace day_bins = 100 if int_days > 120
    
    * Only HD sections
    keep if dep == 1
    
    * Run regressions and store estimates
    * For mixed_highc_hindus
    regress rating c.mixed_highc_hindus#i.day_bins i.line_section mean_contact_hindu_ls i.Shift i.doy
    estimates store highc
    
    * For mixed_lowc_hindus
    regress rating c.mixed_lowc_hindus#i.day_bins i.line_section mean_contact_hindu_ls i.Shift i.doy
    estimates store lowc
    
    * Use suest to combine estimates
    suest highc lowc
    * Note: Cannot use cluster() in the initial regressions when using suest
    * Cluster at the line_section_team level in suest
    suest highc lowc, vce(cluster line_section_team)
    
    * Use parmest to create a new dataset with the estimates
    parmest, norestore level(90 95)
    
    * Save the dataset
    tempfile estimates
    save `estimates'
    
    // Load the estimates and reshape the data for easier plotting
use `estimates', clear
keep if regexm(parm, "day_bins#c.mixed")
split parm, parse(".")
drop parm
replace parm1 = "0" if parm1 == "0b"
destring parm1, gen(day_bins)
replace day_bins = day_bins + 2 if parm3 == "mixed_highc_hindus"

// Create the plot
twoway (rspike min90 max90 day_bins if parm3 == "mixed_highc_hindus", msymbol(circle) color(orange%90) lwidth(medthick)) ///
       (rspike min90 max90 day_bins if parm3 == "mixed_lowc_hindus", msymbol(square) color(navy%90) lwidth(medthick)) ///
	   (rspike min95 max95 day_bins if parm3 == "mixed_highc_hindus", msymbol(circle) color(orange%70) lwidth(thin)) ///
       (rspike min95 max95 day_bins if parm3 == "mixed_lowc_hindus", msymbol(square) color(navy%70) lwidth(thin)) ///
       (scatter estimate day_bins if parm3 == "mixed_highc_hindus", msymbol(circle) color(orange)) ///
       (scatter estimate day_bins if parm3 == "mixed_lowc_hindus", msymbol(square) color(navy)) ///
       , legend(label(1 "Hindus with > median Muslim teammates") label(2 "Hindus with < median Muslims teammates") ring(0) pos(5) order(1 2)) ///
       ytitle("Effect on section ratings relative to homogeneous teams") xtitle(Day Bins) ///
	                       ylabel(-0.25(0.1)0.25, nogrid) ///
					xlabel(0(20)100, nogrid) ///
					       xlabel(0 "0-20" 20 "21-40" 40 "41-60" 60 "61-80" 80 "81-100" 100 "> 100", nogrid) ///
title("Heteregenous effects by baseline contact of Hindus with Muslims") yline(0)


graph export "$Output/Figures/hetero_HinduExposure_BW.pdf", replace
    
}




*--------------------------------------*
*       Shift Heterogeneity Analysis   *
*--------------------------------------*
if `shift_het' == 1 {

 
    * Load Cleaned Data
    use "$Data/Final/linesectionlevel_datamerged.dta", clear
    
    * Main Regressions (HD and LD pooled)
    	
    * Day Shifts
    reghdfe rating mixedXLD mixedXHD if inlist(Shift, 1, 2), absorb(doy line_section Shift) cluster(line_section_team)
    eststo A1
    local obs1 = e(N)
    summarize rating if e(sample)
    local meanA = string(round(r(mean), .01))
	test mixedXHD = mixedXLD
    
    * Night Shift
    reghdfe rating mixedXLD mixedXHD if inlist(Shift, 3), absorb(doy line_section Shift) cluster(line_section_team)
    eststo A2
    local obs2 = e(N)
    summarize rating if e(sample)
    local meanB = string(round(r(mean), .01))
    test mixedXHD = mixedXLD

   
    *--------------------------------------*
    *             Output Table             *
    *--------------------------------------*
    
    estout A1 A2 using "$Output/Tables/tables_rating_HDLD_shifthet", style(tex) replace ///
        keep(mixedXLD mixedXHD) ///
        cells(b(star fmt(%9.4f)) se(par)) ///
        nolabel collabels(none) mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) ///
        varlabels(mixedXLD "Mixed X LD" ///
                  mixedXHD "Mixed × HD")
    
    * Notes for the table
    local tex " \\ \hline"
    local tex "`tex' Observations & `obs1' & `obs2' \\"
    local tex "`tex' Mean Dep. Var & `meanA' & `meanB' \\"
    local tex "`tex' Day Effects & Yes & Yes  \\"
    local tex "`tex' Shift Effects & Yes & Yes  \\"
    local tex "`tex' Line × Section Effects & Yes & Yes \\"
    local tex "`tex' \multicolumn{5}{p{12cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level. "
    local tex "Mixed is a dummy variable coded 1 if the line-section-level team is religiously mixed. "
    local tex "*** p$<$0.01, ** p$<$0.05, * p$<$0.1.} \\ \end{tabular} }"
    
    esttab A1 A2  using "$Output/Tables/tables_rating_HDLD_shiftheta", style(tex) replace booktabs ///
        d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
        mtitles("(1)" "(2)") ///
        mgroups("Day Time Shifts" "Night Time Shifts", pattern(1 1) span) ///
        title("Shift Heterogeneity Analysis") ///
        alignment(D{.}{.}{-1})
		
}





*--------------------------------------*
* Heterogeneous Effects by Baseline    *
* Political Attitudes of Hindus        *
*--------------------------------------*

if `bjp_nrc_hindu' == 1 {

    *--------------------------------------*
    *         Set Up Environment           *
    *--------------------------------------*

    set more off
    clear all
    set segmentsize 3g
    set scheme plotplainblind


    * Load Cleaned Data
    use "$Data/Final/linesectionlevel_datamerged.dta", clear

    *--------------------------------------*
    * Create Variables Based on Political  *
    *        Attitudes of Hindus           *
    *--------------------------------------*

    * Calculate the median of political attitudes among Hindus
    quietly summarize pattitudes_hindu if !missing(pattitudes_hindu), detail
    local median_attitude = r(p50)

    * Create variables for high and low Hindutva support among Hindus in mixed teams
    generate mixed_high_support = mixed
    replace mixed_high_support = . if pattitudes_hindu <= `median_attitude' & mixed == 1

    generate mixed_low_support = mixed
    replace mixed_low_support = . if pattitudes_hindu > `median_attitude' & mixed == 1

    *--------------------------------------*
    *            Create Day Bins           *
    *--------------------------------------*

    egen day_bins = cut(int_days), at(0(20)120)
    replace day_bins = 100 if int_days > 120

    * Only High Dependency (HD) sections
    keep if dep == 1

    *--------------------------------------*
    *      Run Regressions and Store       *
    *             Estimates                *
    *--------------------------------------*

    * Regression for mixed_high_support (High Hindutva Support)
    regress rating c.mixed_high_support#i.day_bins i.line_section i.Shift i.doy
    estimates store high_support

    * Regression for mixed_low_support (Low Hindutva Support)
    regress rating c.mixed_low_support#i.day_bins i.line_section i.Shift i.doy
    estimates store low_support

    *--------------------------------------*
    *    Combine Estimates and Prepare     *
    *             for Plotting             *
    *--------------------------------------*

    * Use suest to combine estimates and specify clustering
    suest high_support low_support, cluster(line_section_team)

    * Use parmest to extract parameters for plotting
    parmest, norestore level(90 95)

    * Save the estimates
    tempfile estimates
    save `estimates'



    *--------------------------------------*
    *             Create Plot              *
    *--------------------------------------*

    // Load the estimates and reshape the data for easier plotting
    use `estimates', clear
    keep if regexm(parm, "day_bins#c.mixed")
    split parm, parse(".")
    drop parm
    replace parm1 = "0" if parm1 == "0b"
    destring parm1, gen(day_bins)
    replace day_bins = day_bins + 2 if parm3 == "mixed_high_support"

     // Create the plot
      twoway (rspike min90 max90 day_bins if parm3 == "mixed_high_support", msymbol(circle)  color(orange%90) lwidth(medthick)) ///
       (rspike min90 max90 day_bins if parm3 == "mixed_low_support", msymbol(square) color(navy%90) lwidth(medthick)) ///
	   (rspike min95 max95 day_bins if parm3 == "mixed_high_support", msymbol(circle)  color(orange%70) lwidth(thin)) ///
       (rspike min95 max95 day_bins if parm3 == "mixed_low_support", msymbol(square) color(navy%70) lwidth(thin)) ///
       (scatter estimate day_bins if parm3 == "mixed_high_support", msymbol(circle) color(orange)) ///
       (scatter estimate day_bins if parm3 == "mixed_low_support", msymbol(square) color(navy)) ///
       , legend(label(1 "Hindus with > median support Hindutva Politics") label(2 "Hindus with < median support Hindutva Politics") ring(0) pos(5) order(1 2)) ///
       ytitle("Effect on section ratings relative to homogeneous teams") xtitle(Day Bins) ///
	                 ylabel(-0.25(0.1)0.25, nogrid) ///
					xlabel(0(20)100, nogrid) ///
					       xlabel(0 "0-20" 20 "21-40" 40 "41-60" 60 "61-80" 80 "81-100" 100 "> 100", nogrid) ///
title("Heteregenous effects by baseline political attitudes of Hindus") yline(0)

graph export "$Output/Figures/hetero_PoliAttitude_BW.pdf", replace

}


*--------------------------------------*
* Heterogeneous Effects by Baseline    *
* Exposure Index                       *
*--------------------------------------*

if `baselineexp_index' == 1 {

    *--------------------------------------*
    *         Set Up Environment           *
    *--------------------------------------*

    set more off
    clear all
    set segmentsize 3g
    set scheme plotplainblind

  

    * Load Cleaned Data
    use "$Data/Final/linesectionlevel_datamerged.dta", clear

    *--------------------------------------*
    *   Create Variables Based on Baseline *
    *           Exposure Index             *
    *--------------------------------------*

    * Calculate the median of baseline exposure index
    quietly summarize baselineexp_index if !missing(baselineexp_index), detail
    local median_exposure = r(p50)

    * Create variables for high and low baseline exposure in mixed teams
    generate mixed_high_exposure = mixed
    replace mixed_high_exposure = . if baselineexp_index <= `median_exposure' & mixed == 1

    generate mixed_low_exposure = mixed
    replace mixed_low_exposure = . if baselineexp_index > `median_exposure' & mixed == 1

    *--------------------------------------*
    *            Create Day Bins           *
    *--------------------------------------*

    egen day_bins = cut(int_days), at(0(20)120)
    replace day_bins = 100 if int_days > 120

    * Only High Dependency (HD) sections
    keep if dep == 1

    *--------------------------------------*
    *      Run Regressions and Store       *
    *             Estimates                *
    *--------------------------------------*

    * Regression for mixed_high_exposure (High Baseline Exposure)
    regress rating c.mixed_high_exposure#i.day_bins i.line_section i.Shift i.doy
    estimates store high_exposure

    * Regression for mixed_low_exposure (Low Baseline Exposure)
    regress rating c.mixed_low_exposure#i.day_bins i.line_section i.Shift i.doy
    estimates store low_exposure

    *--------------------------------------*
    *    Combine Estimates and Prepare     *
    *             for Plotting             *
    *--------------------------------------*

    * Use suest to combine estimates and specify clustering
    suest high_exposure low_exposure, vce(cluster line_section_team)

    * Use parmest to extract parameters for plotting
    parmest, norestore level(90 95)

    * Save the estimates
    tempfile estimates
    save `estimates'

    *--------------------------------------*
    *         Prepare Data for Plot        *
    *--------------------------------------*

    // Load the estimates and reshape the data for easier plotting
use `estimates', clear
keep if regexm(parm, "day_bins#c.mixed")
split parm, parse(".")
drop parm
replace parm1 = "0" if parm1 == "0b"
destring parm1, gen(day_bins)
replace day_bins = day_bins + 2 if parm3 == "mixed_high_exposure"

// Create the plot
twoway (rspike min90 max90 day_bins if parm3 == "mixed_high_exposure", msymbol(circle) color(orange%90) lwidth(medthick)) ///
       (rspike min90 max90 day_bins if parm3 == "mixed_low_exposure", msymbol(square) color(navy%90) lwidth(medthick)) ///
       (rspike min95 max95 day_bins if parm3 == "mixed_high_exposure", msymbol(circle) color(orange%70) lwidth(thin)) ///
       (rspike min95 max95 day_bins if parm3 == "mixed_low_exposure", msymbol(square) color(navy%70) lwidth(thin)) ///
       (scatter estimate day_bins if parm3 == "mixed_high_exposure", msymbol(circle) color(orange)) ///
       (scatter estimate day_bins if parm3 == "mixed_low_exposure", msymbol(square) color(navy)) ///
       , legend(label(1 "High Baseline Exposure") label(2 "Low Baseline Exposure") ring(0) pos(5) order(1 2)) ///
       ytitle("Effect on section ratings relative to homogeneous teams") xtitle(Day Bins) ///
	                       ylabel(-0.25(0.1)0.25, nogrid) ///
					xlabel(0(20)100, nogrid) ///
					       xlabel(0 "0-20" 20 "21-40" 40 "41-60" 60 "61-80" 80 "81-100" 100 "> 100", nogrid) ///
title("Heteregenous effects by general (non-work) exposure: HD sections") yline(0)

graph export "$Output/Figures/hetero_non_work_Exposure_BW.pdf", replace

}


/*---------------------------------*
 *      Spillover Analysis         *
 *---------------------------------*/

if `spillover' == 1{
	
	*Load Data
	 use "$Data/Final/linesectionlevel_datamerged.dta", clear
	
	
	

/* Regression 1: Table B.24 - Model A1 */
areg rating mixed i.line_section i.Shift ///
    if inlist(Line, "Old Layer", "Layer Line"), ///
    absorb(doy) cluster(line_section_team)
eststo A1

/* Capture statistics for Model A1 */
local obs1 = e(N)
quietly summarize rating if e(sample)
local meanA = string(round(r(mean), 0.01))

/* Regression 2: Table B.24 - Model A2 */
areg rating c.mixed#i.dep i.line_section i.Shift ///
    if inlist(Line, "Old Layer", "Layer Line"), ///
    absorb(doy) cluster(line_section_team)
eststo A2

/* Capture statistics for Model A2 */
local obs2 = e(N)
quietly summarize rating if e(sample)
test 0.dependency#c.mixed = 1.dependency#c.mixed
local p1 = string(round(r(p), 0.01))
local meanB = string(round(r(mean), 0.01))

/* Export Table B.24 */
estout A1 A2 using "$Output/Tables/tables_rating_LDHD.tex", style(tex) replace ///
    keep(mixed 0.dependency#c.mixed 1.dependency#c.mixed) ///
    cells(b(star fmt(%9.4f)) se(par)) ///
    nolabel collabels(none) mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) ///
    varlabels(mixed "Mixed" ///
        0.dependency#c.mixed "Mixed X LD" ///
        1.dependency#c.mixed "Mixed X HD")

/* Construct LaTeX footer for Table B.24 */
local tex "\hline"
local tex "`tex' Observations & `obs1' & `obs2' \\"
local tex "`tex' Mean Dep. Var & `meanA' & `meanB' \\"
local tex "`tex' Education and Tenure Controls & Yes & Yes \\"
local tex "`tex' Day Effects & Yes & Yes \\"
local tex "`tex' Shift Effects & Yes & Yes \\"
local tex "`tex' Line*Section Effects & Yes & Yes \\"
local tex "`tex' \multicolumn{5}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level. Mixed is a dummy variable coded 1 if the line-section-level team is religiously mixed. *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"

/* Additional Export with esttab for Table B249 */
esttab A1 A2 using "$Output/Tables/tables_rating_HDLDa.tex", style(tex) replace booktabs ///
    d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
    mtitles("(1)" "(2)") ///
    mgroups("Ratings" "Ratings", pattern(1 1) ///
    prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))

/*---------------------------------*
 *      Spillover Analysis B.23     *
 *---------------------------------*/

areg rating c.mixed#i.dep i.line_section i.Shift if section_n <= 1, absorb(doy) cluster(line_section_team)
 eststo A1
 local obs1 = e(N)
 sum rating if e(sample) == 1
 local meanA = string(round(r(mean), .01))


 areg rating_dum c.mixed#i.dep i.line_section i.Shift if section_n <= 1, absorb(doy) cluster(line_section_team)
 eststo A2
 local obs2 = e(N)
 su rating if e(sample) == 1
 local meanB = string(round(r(mean), .01))
 
 areg rating c.mixed#i.dep i.line_section i.Shift if section_n <= 2, absorb(doy) cluster(line_section_team)
 eststo A3
 local obs3 = e(N)
 test 0.dependency#c.mixed = 1.dependency#c.mixed
 local p1 = string(round(r(p), .01))
 sum rating if e(sample) == 1
 local meanC = string(round(r(mean), .01))
 
 areg rating_dum c.mixed#i.dep i.line_section i.Shift if section_n <= 2, absorb(doy) cluster(line_section_team)
 eststo A4
 local obs4 = e(N)
 test 0.dependency#c.mixed = 1.dependency#c.mixed
 local p2 = string(round(r(p), .01))
 sum rating if e(sample) == 1
 local meanD = string(round(r(mean), .01))


 
 estout A1 A2 A3 A4 using  "$Output/Tables/tables_rating_pooled4", style(tex) replace ///
	keep(0.dependency#c.mixed  1.dependency#c.mixed) ///
	cells(b(star fmt(%9.4f)) se(par)) ///
	nolabel collabels(none) mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) ///
	varlabels( 0.dependency#c.mixed  "Mixed X LD" ///
				1.dependency#c.mixed  "Mixed X HD")
				 
local tex " \\ \hline"
local tex "`tex' Observations & `obs1' & `obs2' & `obs3' & `obs4'  \\"
local tex "`tex'p(Mixed X HD = Mixed X LD) &  &  & `p1' & `p2'  \\"
local tex "`tex' Mean Dep. Var & `meanA' & `meanB' & `meanC' & `meanD'  \\"
local tex "`tex' Day Effects & Yes & Yes & Yes & Yes  \\"
local tex "`tex' Shift Effects & Yes & Yes & Yes & Yes \\"
local tex "`tex' Line*Section Effects & Yes & Yes & Yes & Yes  \\"
local tex "`tex' \multicolumn{5}{p{8cm}}{\tiny \textit{Notes:} Standard errors clustered at line-section-team level."
local tex "`tex' Mixed  is  a  dummy  variable  coded  1  if  the  line-section-level  team  is  religiously mixed." 
local tex "`tex' *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"
	
esttab A1 A2 A3 A4 using "$Output/Tables/tables_rating_pooled4a", style(tex) replace booktabs ///
	d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
	mtitles("(1)" "(2)" "(3)" "(4)") ///
	mgroups("Ratings" "Ratings > Median" "Ratings" "Ratings > Median", pattern(1 1 1 1) ///
	prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
}


/*--------------------------------------------------*
 *      Robustness to Alternate HD/LD Thresholds     *
 *-------------------------------------------------*/

if `dep_robustness' == 1{

* Load  Data
use "$Data/Final/linesectionlevel_datamerged.dta", clear

* Merge with dependency data
merge m:1 Line section using "$Data/Original/line_section_dependency_mins.dta"
drop _merge

gen dep_alternative = w_int_mins > 4 if !mi(w_int_mins)

 areg rating c.mixed#i.dep_alternative i.line_section i.Shift, absorb(doy) cluster(line_section_team)
 eststo A1
 local obs1 = e(N)
 test 0.dep_alternative#c.mixed = 1.dep_alternative#c.mixed
 local p1 = string(round(r(p), .01))
 sum rating if e(sample) == 1
 local meanA = string(round(r(mean), .01))
 drop dep_alternative
 
 gen dep_alternative = w_int_mins > 6 if !mi(w_int_mins)

 areg rating c.mixed#i.dep_alternative i.line_section i.Shift, absorb(doy) cluster(line_section_team)
 eststo A2
 local obs2 = e(N)
 test 0.dep_alternative#c.mixed = 1.dep_alternative#c.mixed
 local p2 = string(round(r(p), .01))
 sum rating if e(sample) == 1
 local meanB = string(round(r(mean), .01))
 drop dep_alternative
 
  gen dep_alternative = w_int_mins > 8 if !mi(w_int_mins)

 areg rating c.mixed#i.dep_alternative i.line_section i.Shift, absorb(doy) cluster(line_section_team)
 eststo A3
 local obs3 = e(N)
 test 0.dep_alternative#c.mixed = 1.dep_alternative#c.mixed
 local p3 = string(round(r(p), .01))
 sum rating if e(sample) == 1
 local meanC = string(round(r(mean), .01))
 drop dep_alternative
 
  estout A1 A2 A3 using  "$Output/Tables/tables_rating_depthresh", style(tex) replace ///
	keep( 0.dep_alternative#c.mixed  1.dep_alternative#c.mixed) ///
	cells(b(star fmt(%9.4f)) se(par)) ///
	nolabel collabels(none) mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) ///
	varlabels( 0.dep_alternative#c.mixed  "Mixed X LD" ///
				1.dep_alternative#c.mixed  "Mixed X HD")
 
 local tex " \\ \hline"
local tex "`tex' Observations & `obs1' & `obs2' & `obs3'   \\"
local tex "`tex'p(Mixed X HD = Mixed X LD) &  `p1' & `p2' &  `p3' \\"
local tex "`tex' Mean Dep. Var & `meanA' & `meanB' & `meanC'  \\"
local tex "`tex' Day Effects & Yes & Yes & Yes \\"
local tex "`tex' Shift Effects & Yes & Yes & Yes  \\"
local tex "`tex' Line*Section Effects & Yes & Yes & Yes  \\"
local tex "`tex' \multicolumn{4}{p{8cm}}{\tiny \textit{}"
local tex "`tex' Mixed  is  a  dummy  variable  coded  1  if  the  line-section-level  team  is  religiously mixed." 
local tex "`tex' *** p<0.01, ** p<0.05, * p<0.1.} \\ \end{tabular} }"
	
esttab A1 A2 A3 using "$Output/Tables/tables_rating_depthresha", style(tex) replace booktabs ///
	d(*) nolabel collabels(none) noobs postfoot("`tex'") nonum ///
	mtitles("(1)" "(2)" "(3)") ///
	mgroups("HD (> 4 mins)" "HD (> 6 mins)"  "HD (> 8 mins)", pattern(1 1 1) ///
	prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}))
	
 
}

